!======================================================
! A. Gazizov, LNGS, INFN, 04.12.10
! Contains common values to pass to functions
!======================================================
! January 2011, LNGS
!======================================================
! A. Gazizov, LNGS, update on 15.01.11

module commirr
   implicit none
!-------------------------------------------------------------------------------------------------
   real(8),parameter:: Omega_m = 0.3d0 & ! cosmology
                     , Lambda  = 0.7d0 & !
                     , ln10    = 2.3025850929940457d0 ! ln(10)
   real(8) :: z_f, z_f1, lnzf1         & ! max redshift of epoch  z_f ans 1+z_f, ln(z_f1) 
            , tmax                     & ! maximum time \tau = ln(z_f1)
            , gam                      & ! \gamma_g
            , evl                      & ! Q*(1+z)^evl 
            , E_0, E_min               & 
            , lgE_max                  & ! E in eV
            , E_max, E_c               &
            , Emax, EE0                &
            , lgE0                     &
            , xmin, xmax, deltax      

!   integer, parameter:: npdesp=2, ldy=npdesp  !, nxp=2001
   integer,parameter:: nxp    = 5001   &     
                     , npdesp = 2      &
                     , ldy    = npdesp &
                     , nstep  = 61
                      
   integer, private:: i, j
   real(8):: Xmir   & ! T_mirr/T
           , Xmir3    ! Xmirr^3

! Acording to ZB 13.01.11
! R_E = 0.5/(1 + \omega \Gamma_n \tau_{\rm osc})  
! \omega = 4.5\cdot 10^{-11}\left(\frac{B}{10^{-14} G} \right) sec^{-1} 
! \tau_{\rm osc} = 1 sec

   real(8):: B_15, tau_osc
   real(8):: r_crit  ! ln(E_crit/E_0)    
   real(8):: r_0
   real(8),dimension(nxp    ):: xbreak, dlgEeV
   real(8),dimension(ldy,nxp):: dlgE3J, y
   real(8),dimension(nstep  ):: tgrid

   integer:: ido 
   real(8):: tol, hinit   ! parameters for DMOLCH 

   character( 1):: Osc
   character( 2):: chzf
   character( 3):: chgam
   character( 4):: chEmax
   character( 4):: chevl, chXmir, chB15
   character(28):: fname

   character( 1):: ohne

   integer:: jbegin    ! to resume calculations

!-------------------------------------------------------------------------------------------------
                        CONTAINS
!-----------------------------------------------------------------------------

 Subroutine ReadParamMirr
   implicit none

   open(1,file='ParamMirr.dat',status='old')
      read(1,'(10X,d7.1)') tol
      read(1,'(10X,d7.1)') E_0
      read(1,'(10X,d7.1)') E_min
      read(1,'(10X,d6.1)') tau_osc
      read(1,'(11X,A1  )') Osc
      read(1,'(10X,d6.1)') z_f
      read(1,'(11X,A2  )') chzf
      read(1,'(10X,d6.1)') gam
      read(1,'(11X,A3  )') chgam
      read(1,'(10X,d7.2)') lgE_max
      read(1,'(11X,A4  )') chEmax
      read(1,'(10X,d6.1)') evl
      read(1,'(11X,A4  )') chevl
      read(1,'(10X,d6.1)') Xmir
      read(1,'(11X,A4  )') chXmir
      read(1,'(10X,d6.1)') B_15
      read(1,'(11X,A4  )') chB15
      read(1,'(10X,I3  )') jbegin
      read(1,'(11X,A1  )') ohne

   print*, ' File ParamMirr.dat is read!'
 
   if(Osc == 'N') then
     fname = Osc//chzf//chgam//chEmax//chevl//'.dat'
   else
     fname = Osc//chzf//chgam//chEmax//chevl//chXmir//'.dat'
   endif
   print*, ' '//fname

   E_max = 1d1**lgE_max

!===
   hinit = 1d-2*tol   ! initial step
   ido   = 1          ! initial value of ido
!===

   z_f1  = 1d0 + z_f  ! 1 + z_f
   lnzf1 = dlog(z_f1) ! ln(1 + z_f)
   tmax  = lnzf1
   do j = 1, nstep
      tgrid(j) = tmax*(dble(j-1)/dble(nstep-1))**2
   enddo

   if(B_15 == 0d0) then
      r_crit = 1d4
   else
      r_crit = dlog(2d20/tau_osc/B_15/E_0)
   endif
   
!  The range of energies
   r_0 = dlog(E_0)          ! ln(E_0/eV)
   E_c = E_0/1d9            ! in GeV
   Emax = E_max/1d9         ! in GeV 

! Recalculated range of x=ln(E/E_0)
   xmin = dlog(E_min/E_0)   ! ln(E_min/E_0)	 (ln( 0.1) for E_min = 10^{17} eV )
   xmax = dlog(E_max/E_0)   ! ln(E_max/E_0)  (ln(1000) for E_max = 10^{22} eV )

   deltax = xmax - xmin

   Xmir3 = Xmir*Xmir*Xmir

!  xbreak = ln(E/E_0),  E_0 = 10^18 eV         
   
   do i = 1,nxp
      xbreak(i) = xmin + deltax*dble(i-1)/dble(nxp-1)
      dlgEeV(i) = (xbreak(i) + r_0)/ln10  ! lg(E/eV)
   enddo

   y = 0d0

 end Subroutine ReadParamMirr
!-----------------------------------------------------------------------------

!***************************************************
! Linear interpoalation

   real(8) Function FLin(x,xa,xb,ya,yb)
      implicit none
      real(8),intent(in):: x, xa, xb, ya, yb

      if(xa == xb) then
         print*, ' Stop, xa = xb in FLin !!!'
         stop
      else
         FLin = ya + (yb - ya)/(xb-xa)*(x-xa)
      endif   
   end Function FLin
!***************************************************	 


end module commirr
!==============================================================


